A mathematical model and numerical simulation for SARS-CoV-2 dynamics

Since its outbreak the corona virus-19 disease has been particularly aggressive for the lower respiratory tract, and lungs in particular. The dynamics of the abnormal immune response leading to lung damage with fatal outcomes is not yet fully understood. We present a mathematical model describing the dynamics of corona virus disease-19 starting from virus seeding inside the human respiratory tract, taking into account its interaction with the components of the innate immune system as classically and alternatively activated macrophages, interleukin-6 and -10. The numerical simulations have been performed for two different parameter values related to the pro-inflammatory interleukin, searching for a correlation among components dynamics during the early stage of infection, in particular pro- and anti-inflammatory polarizations of the immune response. We found that in the initial stage of infection the immune machinery is unable to stop or weaken the virus progression. Also an abnormal anti-inflammatory interleukin response is predicted, induced by the disease progression and clinically associated to tissue damages. The numerical results well reproduce experimental results found in literature.

. (6) . (7) The above system of PDEs has been integrated using the FEM and imposing zero-flux boundary conditions on the domain boundaries, assuming that all the components stay confined within the simulated biological domain. The two-dimensional spatial domain of 1 mm 2 area has been discretized with a uniform mapped mesh with 2500 domain elements. The time discretization has been done using the implicit Euler method, and the local variable distribution was interpolated with quadratic shape functions via the Galerkin's method, see for example Amoddeo 1,2 , while the FEM is a numerical technique that has been found to be efficient to solve systems of coupled non-linear PDE resulting from complex systems modelling such as fast transitions in nematic liquid crystals 3 .
The PDE system has been non-dimensionalized introducing a set of reference quantities representative of typical biological values 2 , used to rescale the variables and the parameters, as detailed in the following.

Parameters and reference quantities
The parameters were harvested from the literature, when possible, or calibrated so that the simulation results would match experimental results found in the literature. Supplementary Table   S1 is divided into two sections. In the first one we report the reference quantity used for the parameters non-dimensionalization, their symbols, units, values and sources. Multiplying, respectively, x, or y, by l, and t by t, it can recover the spatial and temporal dimensional scales, hence t×D = l 2 .
In the second section we report the parameters used: for each one, starting from the first column, we insert a short description, the used symbol, units, how to obtain the non-dimensional value, the nondimensional value and finally the reference source, when possible. Supplementary Table S1. Summary of the reference quantities and parameters used in the model.
In the following we detail the parameters values or, in the absence, their estimation/assumption or assimilation to other parameters.

Reference Quantities
Reference quantity for viral load -Vr. Sender et al. 4 have estimated the number of virions in humans extrapolating the data collected from experiments on rhesus macaques 2-4 days after SARS-CoV-2 inoculation. In lungs, they estimate a number of 10 6 -10 8 n g -1 . Using then a lung density r=1g/1mL, we extrapolate a range 10 6 -10 8 n cm -3 , while we assume Vr = 10 7 n cm -3 .
Reference quantity for CD4 + -T cells -Tr. It is commonly recognized that in adults and adolescents the number of such cells in the blood can vary in the range 500-1200 cell mm -3 ; we then assume as reference the value 5×10 6 cell cm -3 .
CD4 + -T cells maximum carrying capacity -TM. Set equal to Tr.
Reference quantity for infected CD4 + -T cells -Ir. It is considered the same as Tr.
Reference quantity for M1 activated macrophages -M1r. Leonard et al. 5 assumed the number of macrophages recruited by a tumour lesion is about 6.9×10 6 cell cm -3 . Hence, in the absence of measured data for SARS-CoV-2 infection, we assume the same value.
Reference quantity for M2 activated macrophages -M2r. It is considered the same as M1r.
Reference quantity for IL-6 -Lr. The paper of Liao et al. 6 reports on the mathematical modelling of the anti-tumour response induced by interleukin-27 (IL-27), involving also CD8 + -T cells, IL-10 and IFN-g, where as reference quantity for IL-27, IL-10 and IFN-g it is assumed the value of 10 2 pg cm -3 . Due to its pro-inflammatory role, we assimilate IL-6 to IFN-g, assuming for what concerns us a reference quantity for IL-6 = 10 2 pg cm -3 . Being the molecular mass of IL-6, MIL-6 ≈ 21 kDa 7 , the Avogadro number NA = 6.022×10 23 , then the number of IL-6 cells contained in 10 2 pg is given by ." × 6.022 × 10 24 ≈ 2.87 × 10 5 .
Reference quantity for IL-10 -Nr. Still from the previous case, a reference quantity for IL-10 it is assumed 6 as 10 2 pg cm -3 . Being the molecular mass of IL-10, MIL-10 ≈ 17-21 kDa 8 , choosing the higher value for convenience, as for Lr we obtain Nr = 2.87 × 10 5 cell cm -3 .
Characteristic length -l, characteristic diffusion coefficient -D, characteristic time scale -t.
Our simulations are modelled in a square domain with side 1 mm, which is a typical length for biological phenomena, then l = 0.1 cm.
CD4 + -T cell diffusion coefficient -DT. In Lai et al. 10  CD4 + -T cells infection rate -k. In Quirouette et al. 9 an estimate for the infection rate of target cell by influenza A virus is given as 1.33 × 10 -6 (TCID50 mL -1 ) -1 h -1 . In order to find a conversion from CD4 + -T cells activation rate by IL-6 -f21, T cell inhibition rate by IL-10 -f23. In the absence of direct measures, we adopt by assimilation some parameter used in Liao et al. 6 for tumour cells, having regard to their size and magnitude. In particular, the death rate of tumour cells and the inhibition rate for tumour cells from IL-10, respectively 1.73 × 10 -1 day -1 and 3.45 × 10 -1 day -1 , give values spanning in the 2 × 10 -6 s -1 ÷ 4 × 10 -6 s -1 range. Therefore, we set f21 = 2 × 10 -6 s -1 and f23 = 4 × 10 -6 s -1 .
CD4 + -T cells production rate by M1 -f22, T cell inhibition rate by M2 -f24. These parameters where calibrated in order to match as close as possible the experimental findings in Huang et al. 14 , Zhou et al. 15 and Du et al. 16 , then we set f22 = 1.67 × 10 3 s -1 and f24 = 6.9 s -1 .
Infected CD4 + -T cells decay rate -d. We found in literature 13 that a decline rate for infected cells producing virus for HIV infected T cells is 0.39 day -1 = 4.5 × 10 -6 s -1 , while in Lai et al. 10 a nondimensional value of 0.2 for the same parameter is given. It follows that the decay rate for infected T cells can span in the 10 -6 s -1 ÷ 10 -5 s -1 range: we set d = 2 × 10 -6 s -1 .
Infected CD4 + -T cells reduction rate due to hyper-inflammation -f32. This parameter has been calibrated in order to match as close as possible the experimental findings in Longhi et al. 17 , Huang et al. 14 , Zhou et al. 15 and Du et al. 16 , then we set f32 = 3.5 × 10 -13 cm 3 cell -1 s -1 .
IL-10 inhibition rate induced by IL-6 -f73. Nickaeen et al. 18 proposed an agent-based mathematical model for macrophage polarization, where the production/inhibition rates of several cytokines entering in their model have been estimated or harvested from literature. In particular, IL-12 in that model inhibits IL-4 production at a rate 9.26 × 10 -5 , then by assimilating IL-12 to IL-6 and IL-4 to IL-10 we adopt such non-dimensional value which in dimensional form reads as f73 = 9.26 × 10 -9 s -1 .
Nevertheless, it can be seen that all the parameters related to cytokine production/inhibition belong to the 9.26 × 10 -5 ÷ 2 × 10 -3 non-dimensional range, even if it comes in response to activation of several pathways that may or may not be involved in our case, but we trust that the orders of magnitude of the parameters are adaptable to our purposes. Similar estimates have been done in Frieboes et al. 19 , for protein production rate and their ranges of variability, as well for protein decay rate, based on proteomic analysis, in a study aimed at the modelling of production, transport and shedding of proteins in a vascularized tumour. Hence, the remaining parameter estimate, is given accordingly, in order to match above mentioned experimental results 14-17 : M1 macrophages production rate induced by IL-6 -f41.  Table 1. to f52 = 1 × 10 -3 and f52 = 2 × 10 -3 , respectively. All the other parameters are as in Table 1. to f52 = 1 × 10 -3 and f52 = 2 × 10 -3 , respectively. All the other parameters are as in Table 1. days), t = 50 (~5.79 days) and t = 60 (~6.94 days), while the first and second column of panels refer to f52 = 1 × 10 -3 and f52 = 2 × 10 -3 , respectively. All the other parameters are as in Table 1 3.5×10  Table 1.